On the emergence of helicity in rotating stratified turbulence 
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We perform numerical simulations of decaying rotating stratified turbulence and show, in the 
Boussinesq framework, that helicity (velocity- vorticity correlation), as observed in super-cell storms 
and hurricanes, is spontaneously created due to an interplay between buoyancy and rotation common 
to large-scale atmospheric and oceanic flows. Helicity emerges from the joint action of eddies and 
of inertia-gravity waves (with inertia and gravity with respective associated frequencies / and A'^), 
and it occurs when the waves are sufficiently strong. For N/ / < 3 the amount of helicity produced 
is correctly predicted by a quasi-linear balance equation. Outside this regime, and up to the highest 
Reynolds number obtained in this study, namely Re « 10000, helicity production is found to be 
persistent for N/f as large as ~ 17, and for ReFr^ and ReRo^ respectively as large as ~ 100 and 
f» 24000. 

PACS numbers: 47.55.Hd, 47.32. Ef, 47.27.-i, 47.27.ek 



I. INTRODUCTION 



Symmetry breaking is a fundamental concept which 
has been quite fruitful in many physical applications 
For a fluid the simplest way to break symmetry is to 
introduce helicity. In that case, the velocity covariancc 
tensor is still expressible in terms of the magnitude of the 
distance between points, i.e., the fluid still has isotropic 
statistics, but mirror symmetry is broken whereby the co- 
variance matrix has an anti-symmetric component which 
can be shown to be proportional to the total helicity 
Hv = (u • u;), which is defined as the correlation between 
the velocity u and its curl, the vorticity lj. Helical struc- 
tures abound in nature, from macroscopic organisms to 
elastomers; helical structures can cause erosion in river 
bends and alter nutrient mixing properties in estu- 
aries, in particular when interacting with tidal flows Q. 
Helical flows are observed as well in the atmosphere, in 
dust devils, tornadoes and hurricanes [J- 

Helicity is an invariant of motion of the non-dissipative 
fluid equations involving the topology of field lines [5| , in- 
cluding in the presence of rotation, but not with stratifi- 
cation, where instead potential vorticity is invariant and 
is essential in determining structures such as sharp jets in 
planetary atmospheres [6| . Invariants are known to play a 
fundamental role in turbulence, since the nonlinear terms 
have to preserve such invariants at the level of triadic in- 
teractions in the incompressible case. However, helical 
(corkscrew) motions do not seem to alter the dynamics 
of homogeneous isotropic turbulence in the absence of 
both rotation and stratification, with the kinetic energy 
and helicity spectral densities (respectively Ev{k) and 
Hv{k), with / Ev{k)dk = \ (u^) and / Hv{k)dk = Hy) 
both following a Kolmogorov spectrum. This implies a 



slow ^ 1/k decay of the relative helicity in Fourier space 
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with (Ty(x) = cos(u, the degree of alignment between 
velocity and vorticity in configuration space. However, it 
is straightforward to show that helicity is created point- 
wise by the alignment of vorticity and pressure or shear 
gradients 0, and it is observed to be strong (cry ^ ±1) 
in the vortex filaments that are ubiquitous in isotropic 
fiuid turbulence at small scale. 

Invariants are also the stepping stone to determine in- 
ertial range behavior in turbulent flows; this principle 
is at the basis of statistical mechanics that has proven 
useful in predicting, for example, the inverse cascade of 
energy for a two-dimensional fluid Q. In fact, a recent di- 
rect numerical simulation of the ideal three-dimensional 
fluid equations in the absence of waves, showed that, 
at intermediate times, a Kolmogorov spectrum devel- 
ops at large scale, the effective dissipation for the large- 
scale fluid being produced by the eddy viscosity emanat- 
ing from the small-scale equilibrated modes How- 
ever, non-conserved quantities can also play an impor- 
tant role through other mechanisms such as interactions 
with waves and large-scale hydrostatic and geostrophic 
balance [lo| . 

When the fluid is conducting, magnetic helicity is an 
invariant in the ideal case and is central to minimum 
energy equilibria in plasmas such as in spheromaks, or 
in solar coronal mass ejections It is also known 

that the generation of large-scale magnetic fields occurs 
due to small-scale mechanic helicity Hy , and that in the 
presence of both rotation and stratification, helicity is 
created and thus a dynamo is facilitated in a wide variety 
of astrophysical settings ■ In the context of this work, 
it is important to note that although in this case the 
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mechanic helicity is not an invariant any longer, it still 
plays an essential role in determining the scaling of the 
fields at large scales. 

Rotating stratified turbulence is important in the at- 
mosphere and oceans, playing a crucial role in their dy- 
namics. In the presence of waves, advective nonlinear 
interactions responsible for the complexity of turbulent 
flows have to compete with the waves and an equilib- 
rium can be reached at some scale and broken at oth- 
ers, the best known example perhaps being the differ- 
ence between the Garret-Munk and the Philligs spectra 
in internal waves in lakes or the oceans [l^ [3] , where 
wave coupling in resonant triads leads to mixing (like in 
coastal currents [Il])j_to vertical dispersion jl6j and to 
enhanced dissipation . A particular set of modes plays 
a major role, in the so-called slow manifold for which the 
frequency of the waves tends to zero, and only turbu- 
lent interactions and standing potential vorticity modes 
remain. When rotation (only) is present in the fluid, 
strong relative helicity can alter the scaling of the distri- 
bution of energy among scales and lead to the occurrence 
of helical long-lived structures [l^ ■ 

What happens when stratiflcation is also included? In 
the remainder of this paper, we address the question of 
rotating stratified turbulence in the absence of forcing, 
as studied for example in [l^, |2^ , but concentrating on 
the creation of helicity and on the link between the evo- 
lution of helicity and the balance of forces such as ro- 
tation, stratification, inertia, and pressure gradients, at 
scales large enough that the effect of nonlinearities (iner- 
tia) is small for strong waves. In the limit of zero nonlin- 
earities, the resulting geostrophic balance is crucial for 
weather forecasting and simulations of climate change. 
However, the consequences of the interplay between ro- 
tation and stratification, as far as helical motions are 
concerned, have been mostly ignored except for the pio- 
neering work of Hide 2l[. In spite of this, helicity was 



hypothesized to be important in the atmosphere in the 
dynamics and persistence of rotating convective storms 
[22l | on the basis of the weakening of non-linear interac- 
tions in the so-called Lamb vector ux uj. 

It is also interesting that helicity is measured in the 
context of forecasting storms and tropical tornadoes, in 
particular in the presence of strong shear, and it can be 
used as an indicator of storm occurrence [23| . Note that 
it has been shown that in some cases (using a specific 
fully helical Beltrami forcing function), shear is created 
at large-scale in a rotating flow [IJl . 

Since helicity in rotating and stratifled flows is no 
longer an invariant even in the absence of dissipation, its 
presence in these atmospheric storms can be accounted 
for but the physical mechanisms governing its creation, 
and the structures associated with it, remain unclear. In 
this paper, we perform a parametric study using direct 
numerical simulations in which we vary both rotation and 
stratification. In that framework, we show that a strongly 
rotating stratifled flow can spontaneously create helicity 
at large scales. 



II. EQUATIONS AND NUMERICAL 
PROCEDURE 



A. Boussinesq equations and parameters 

We integrate the incompressible Boussinesq equations 
in the rotating frame, with constant (solid body) rotation 
r2 and gravity g, anti-aligned in the vertical (z) direction, 
with 9 the buoyancy (in units of velocity), w the verti- 
cal velocity, P the pressure, v the viscosity, and k the 
diffusivity: 

dtU + u Vu- i/Au = - VP - N9ez - 2Qe^ x u, (1) 
dt9 + u-y0- kAO = Nw, (2) 
V • u = . (3) 

We write u = {u, v, w) and we take a unit Prandtl 
number, v ~ k. The Brunt- Vaisala frequency is iV = 
[—gdzO/9Y^^ where 9 is the background imposed stratifi- 
cation. In the general case, one has inertia-gravity waves 
of frequency 



with f :^2n (see, e.g., Bllg). 

The Froudc, Rossby and Reynolds numbers are de- 
fined, respectively, as 

^ ""On "0 „ UqLq 

Fr ^ — — , Ro = — - , Re = , 

NLq fLo V 

with ?io = 1 and Lq = 27r/fco respectively the r.m.s. ve- 
locity and the scale of the initial conditions. These pa- 
rameters vary in the range 0.0063 < Fr < 0.2, 0.0063 < 
Ro < 3.24, and Re w 4000 for grids of 256^ points, while 
Re w 10000 using 512'^ points. Decay is left to occur 
for 3.6 to 7.2 turn-over times, tjvl = Lq/uq. The initial 
velocity field is random, with all three components non- 
zero, and it is centered around wavenumbers ko = [1,2]. 
At < = 0, 6* = 0, and Hy ~ +0.2. Other initial values 
have been used as well to ascertain that the results are 
insensitive to them. Note that we do not attempt to take 
initially a balanced flow; for the time-stepping point of 
view, there is no need to do so, since the resolutions we 
employ are high enough that the relatively small Froude 
and Rossby numbers we simulate can be handled with 
an explicit time stepping resolving the smallest eddy- 
turnover time and the smallest Brunt- Vaisala, inertial, 
and inertia-gravity frequencies. Furthermore, the gener- 
ation of gravity waves that compete with turbulent ed- 
dies is part of the overall dynamics of such flows as the 
Reynolds number increases. 

In the ideal (i/ = 0) case, potential vorticity 



PV 



-fN + fd;,9 - Nuj,, +LJ-V9 



is a point- wise invariant, and the total (kinetic plus po- 
tential) energy Et = Ey + Ep is conserved as well (with 
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FIG. 1: E{k)/k^, where E{k) is the kinetic energy spectrum, averaged for one eddy turn-over time after the peak of enstrophy, 
and where p is chosen to compensate for either a wave turbulence law {p = —2, dashed line) or a Kolmogorov law {p = —5/3, 
solid line); k is the isotropic wavenumber. Runs were computed on grids of 512"^ points with initial conditions at ko £ [1,2], 
and Reynolds numbers Re ~ 10000. The dash-dotted line indicates the best fit in the inertial range. Left: N/f — 2.99, 7?, -r. = 
ReRo'^ « 3.84, TZb = ReFr'^ ^ 0.43, with Fr ^ 0.0063 and Ro ^ 0.019. Right: N/f = 4.0, ReRo'^ 1749, ReFr'^ ^ 109, with 
Fr ~ 0.1 and Ro ~ 0.4. Note the steeper spectrum for moderate N/f and low TZb,tz, and a scaling close to a Kolmogorov law 
for larger N/f and substantially larger TZb,tz- 
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FIG. 2: Top and middle: Temporal evolution of the total 
helicity Hv in several runs with different values of Fr, Ro, 
and A'^, as given by the labels. Note that the time averaged 
value of Hv is negative, indicating negative helicity prevails 
in these runs even when the initial value of the helicity is 
positive. On top are runs with the same Fr whereas in the 
middle, runs with N/f = 1 but with different Fr are shown. 
Bottom: Time evolution of the kinetic enstrophy Zv in runs 
with Fr « 0.01 and A'^ = 12.56, and with different values of 
Ro. In all panels, oscillations are due to gravity waves, with 
their period proportional to A''. 



Ep = i {&^)), with respective enstrophies (proportional 
to dissipation when i> ^ and k ^ 0), 

Zv = {io^) , Zp = {\ve\^). 

Note that PV is quadratic and thus its L2 norm is not 
conserved in general by the truncation of Fourier space 
used in any spectral method; however, the nonlinear term 
u) ■ \79 can be neglected in the presence of strong rotation 
and stratification [2BI, resulting in a quantity whose L2 
norm is conserved after truncation of Fourier space. 



B. The GHOST code and the runs 

The numerical simulations have been carried out us- 
ing the Geophysical High-Order Suite for Turbulence 
(GHOST) code. GHOST is a pseudo-spectral frame- 
work that hosts a variety of partial differential equation 
(PDF) solvers optimized for studying turbulence in a 
[0, 2tt]^ tri-periodic box, and with 2""^ or 4*''-order ex- 
plicit Runge-Kutta time stepping schemes. Using a cu- 
bic box and an explicit time stepping method allows in 
principle, given the parameters are right, for resolving 
all scales including the Ozmidov scale and beyond, when 
isotropy recovers (see, e.g., [l^l for the purely rotating 
case). A classical 2/3 de-aliasing rule is used, meaning 
that for a given resolution of Up points per dimension, the 
maximum available wavenumber is kmax = np/3. The 
code uses a hybrid MPI/OpenMP parallelization scheme 
[28| (MPI is the Message Passing Interface library, and 
OpenMP stands for Open Multiprocessing, an interface 
to program shared memory environments) . The code also 
has a third level of parallelization with the recent ad- 
dition of support for Graphic Processing Units (GPUs) 



FIG. 3: (Color online) Visualization of the buoyancy 6 in runs with 512^ grids, for Re ^ 10000, Fr = 0.1, and Ro = 0.4 (left) 
and for the same Re, Fr = 0.025, and Ro = 0.05 (right). The vertical direction is indicated by the blue arrow; dark (blue) and 
light (green) strata represent respectively positive and negative variations in 6 around its mean, with sizable fluctuations and 
structuring, and with more turbulent eddies at higher Froude number. 



and accelerators for the Fast Fourier Transforms (FFTs) . 
Note that the MPI communication required to complete 
the multidimensional Fourier transforms is all-to-all. The 
code uses a "slab" (ID) domain decomposition among 
MPI tasks, and OpenMP provides a second level of par- 
allelization within each slab or MPI task. The code can 
compute in double or single precision based on resolu- 
tion. GHOST performance has been tested on a variety 
of platforms, and has been shown to scale linearly up to 
98304 processors, with grids up to 6144'^ points. Data is 
stored at regular intervals and post-processed, both for 
quantitative analysis and visualization, the latter being 
performed using the VAPOR visualization software j29j . 

In Table H] we give the major parameters of the simula- 
tions used in this paper. Note that we have restricted our 
analysis to moderate values of iV//, in particular we have 
for all cases N/ f > 1/2. This is because, in the purely 
rotating case (A^ — !• 0), helicity is exactly conserved and 
thus as one goes into that parameter regime, the cre- 
ation of helicity has to become negligible with decreas- 
ing N at fixed /; furthermore, many geophysical flows 
are dominated by gravity waves except at the largest 
scales. There are studies that show for example that, 
for purely rotating flows, a turbulence regime affected by 
waves develops for Ro < 0.2, whereas at intermediate 
Rossby numbers nonlinear transfer is reduced but the in- 
verse cascade characteristic of the bi-dimensionalization 
of the flow does not take place [s^l. Also, for strong waves 
(strong rotation or stratification), turbulence barely de- 
velops resulting in steep spectra; this is related to the 
value achieved by the so-called buoyancy Reynolds num- 



ber TZis defined below, and the equivalent concept for 
rotating fiows, TZ-jz. Considering this region of interest in 
parameter space, and given the constraints of computing 
in three dimensions without resorting to modeling of the 
small-scales, only a limited exploration of the parameters 
is performed. 



III. RESULTS 

A. Generation of helicity for small nonlinearity 

As mentioned in the introduction, helicity is not con- 
served in a rotating and stratified flow, and thus helicity 
can in principle be created by the flow evolution. In 
this section we briefly show how a balance of the forces 
at large scales can result in net helicity of a preferred 
sign in the flow. We start from the primitive Boussinesq 
equations given above and simplify them using several 
hypotheses. Assuming stationarity, weak nonlinearities 
and small dissipation at large scales, it results that the 
equilibrium level of helicity in rotating stratified turbu- 
lence is proportional to N/f and to the correlation be- 
tween buoyancy and vertical shear. A result consistent 



with this behavior was originally obtained by Hide [21 1. 

We start with the momentum equation, Eq. As 
later we will compute time averages, we will assume the 
system is in a steady state and neglect the time deriva- 
tive. We will also consider viscous effects are small, 
and neglect the dissipative term. Computing the ver- 
tical derivative of the remaining terms and taking the 
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Fr 


Ro 


in 




n„ — 256^ 


1 


Re — 4189 


1 


0.0063 


0.0063 




2 


0.0063 


0.0127 




3 


0.0063 


0.0190 


- 


4 


0.0063 


0.0507 


- 


5 


0.0063 


0.0728 




6 


0.0063 


0.1013 




7 


0.0084 


0.0084 




8 


0.0127 


n 0063 




g 


0.0127 


0.0127 




10 


0.0127 


0.0190 




11 


0.0127 


n 0953 




12 


0.0127 


0.0317 




13 


0.0127 


0.0384 




14 


0.0127 


0.0507 


— 


15 


0.0127 


0.1013 


— 


16 


0.0127 


0.1458 


_ 


17 


0.0127 


0.2111 




18 


0.0253 


0.0253 




19 


0.0253 


0.0507 


_ 


20 


0.0253 


0.1013 


- 


21 


0.0253 


0.2026 


- 


22 


0.0253 


0.2913 




23 


0.0253 


0.4054 


- 


24 


0.0507 


0.0507 


- 


25 


0.0507 


0.1267 




26 


0.1013 


0.4224 


- 


27 


0.1013 


0.8444 




28 


0.1013 


1.1515 


* 


29 


0.1013 


1.6888 


* 


30 


0.1266 


0.1266 




31 


0.2026 


0.2026 




32 


0.2026 


0.6079 




33 


0.2026 


0.8106 


- 


34 


0.2026 


1.6888 




35 


0.2026 


2.3268 




36 


0.2026 


3.2428 






Up = 512^ 




Re = 10649 


37 


0.0063 


0.0127 




38 


0.0063 


0.0190 


* 


39 


0.0127 


0.0190 




40 


0.0127 


0.0317 




41 








42 


0.0127 


0.0633 




43 


0.0127 


0.1013 




44 


0.0253 


0.0507 




45 


0.1013 


0.4053 




TABLE ] 


: List of runs 


analyzed in 


this paper with some 



characteristic parameters: run number ub., linear resolution 
Up, Reynolds Re, Froude Fr and Rossby Ro numbers. A star 
in the "in" column indicates points that are in the scatter 
plot with N/ f < 3, and two stars indicate those in the plot 
with 7^B < 20 or Tin < 20 (see Figs, g] and [SJ. 



dot product of the result with u. we get 

u • (u • Vu) = -u • V9;,P - NwO ~2u-d^ {ne^, x u) , 

(4) 

where the velocity field u was written with Cartesian 
components {u,v,w). 

The last term in this equation is 



2u ■ dz (ile^ X u) = —f{udzV — vdzu) = fH± 



(5) 



where H^^ is part of the total helicity density. Indeed, we 
decompose the hehcity as Hy = {H±) + where the 

brackets denote an average, and where H±^ is the helicity 
density associated with uj^. 



i7_L = 16_L • (V X 



(6) 



and H+ is the remainder, H+ = udyW — vdxW + wojz- 
With strong rotation and stratification, H± ^ and 
H± alone essentially determines the total helicity. For ex- 
ample, measurements of {H±)^ ^ (where the subindices 
_L and z indicates the averages are volume averages per- 
formed in the horizontal and vertical directions) found 
in modeling simulations of hurricanes are seen to be two 
orders of magnitude larger than (H^) [3l|. Note also 
that the H± density is proportional to the so-called (cell- 
relative) environmental helicity, when integrated over the 
vertical (see, e.g., j23j). 

From Eqs. (g]) and ([5]), it follows that: 



u ■ VdzP + \i-dz{n- Vu)] . (7) 



When integrated over volume, the second term vanishes 
for an incompressible flow. The third term is cubic in the 
velocity and in a turbulent flow proportional to e{fLz)~^, 
where e is the energy flux, and Lz a characteristic vertical 
scale. For flows with strong rotation and stratification, 
this quantity is expected to be small (see below). As a 
result, after integration and neglecting the third term, 
wc obtain 



N 
1 



de 

oz 



(8) 



This expression was obtained before by Hide [2l|, in 
a slightly different form after integrating by parts and 
assuming periodic boundary conditions 



N 
1 



(9) 



It should be noted that the original derivation in [21| 
assumes the nonlinear term is zero and that the flow is 
in gcostrophic balance. In that case, from hydrostatic 
balance dzW = and helicity in the flow vanishes. Small 
nonlinearities are crucial to ensure that the second-order 
correlator in Eqs. ([5]) or ([5]) is non-zero. 

We thus conclude that, if nonlinearities are small, the 
production of helicity in strongly rotating stratified tur- 
bulence is such that its equilibrium level is directly pro- 
portional to N / f , and results from an interplay between 
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rotation and stratification. In the hmit of f ^ oo (no 
stratification), helicity is exactly conserved. In the limit 
of — !■ oo, stratification dominates and the evolution 
of helicity can only be governed by the nonlinear terms, 
the buoyancy and the dissipation [32[. Indeed, in that 
case dissipation is known to play a role in the overall 
dynamics, e.g., in the changes of potential vorticity once 
gravity waves start to break [33| . Finally, it is interesting 
that N/f scaling has also been advocated, for example, 
in the context of statistical mechanics of non-dissipative 
geophysical flows [sl ]. 



B. Nonlinear effects 

Small nonlinearities and negligible dissipation are just 
the beginning of the story, with these assumptions bro- 
ken when overturning takes place. For example, it 
is known that in three-dimensional turbulence without 
waves, the rate of energy dissipation can be evaluated 
phcnomcnologically as e Uq/Lq, no matter how high 
the Reynolds number; this has been demonstrated us- 
ing highly-resolved direct numerical simulations [ssj up 
to grids of 4096"^ points (for the case of a coupling to a 
magnetic field, in which case Alfven waves are present 
and interact with the flow, see [scj in two dimensions 
(2D), and [13 in 3D). 

Similarly, there is a vast literature which concerns it- 
self with the weak non-linear coupling of waves, e.g., 
through resonant interactions (ssj . through weak turbu- 
lence theory (see [s^ for the rotating case, an d lipl fo r the 
stratified case) , through turbulence closures [30|, |4l| 
and more recently through asymptotic approaches 
lisj . In all cases, when the rotation or stratification is 
not strong enough, and/or when the Reynolds number is 
high enough, a situation described by both the buoyancy 



and what can be called the inertial Reynolds number 



7^g = ReFr^, Tin = ReRo^ , 



(10) 



being large enough, nonlinear couplings will take place 
between eddies and waves, sufficiently so that the energy 
will be transferred to small scales in a self-similar manner. 
Let us note here that in the following (see also Table 
El), 'T^B,'R. are evaluated at the peak of enstrophy, using 
dynamical variables, i.e., based on the so-called integral 
scale L„t = j[Ev{k)/k]/Ev ■ 

In these cases, that we will broadly call wave turbu- 
lence, the scaling laws for either rotating or stratified 
fiows have been deduced both phcnomcnologically and 
analytically in the framework of the aforementioned clo- 
sures and theories, and a continuous power-law spectrum 
is expected, steeper than a Kolmogorov spectrum, be- 
cause of the weakening of interactions in the presence of 
waves. The energy flux e is also typically reduced, in gen- 
eral by a factor that is proportional to the ratio of the 
eddy turnover time to the timescale of the waves. Thus, 
we arc in the presence of turbulence, but not the classical 
Kolmogorov turbulence, rather a wave turbulence regime 



that breaks down at small scales (beyond the Ozmidov 
scale). Two typical cases of the isotropic energy spectra 
are shown in Fig. [TJ We observe well developed kinetic 
energy spectra, rather steep for N/f < 3 {E{k) ~ fc~^), 
whereas the power law is closer to a Kolmogorov law for 
N/f ~ 4. In both cases, the dissipation scale ry (eval- 
uated using a Kolmogorov spectrum) is barely resolved: 
one finds 77 w 0.041 and fc^ « 154 for the case where 
N/f ^ 2.99, and ?/ w 0.044 with fc,, « 144 for the case 
where N/f = 4.0. Note that in the latter case, the tur- 
bulence is stronger and the spectra are not quite as well 
resolved, but the point of this study is not to examine 
Kolmogorov turbulence; in the presence of waves, one 
can simulate higher Reynolds numbers at a given res- 
olution than in the absence of waves, again the set of 
governing parameters being TZb,k rather than Re. 

Note that we do not intend to perform a detailed 
analysis of wave-mode and vortical-mode perpendicular 
and parallel energy and helicity spectra here, but rather 
show that turbulence and helicity develop in the fiows 
we study. There are several examples in the literature 
of such studies at high resolution for the energy (i.e., 
(^'-L, fc||), see for example [46l - l48j for the purely strat- 
ified case, and [131 for rotating turbulence). In our case, 
the choice of isotropic spectra is sufficient to show that 
there is indeed for these parameter regimes power-law 
spectra that develop through nonlinear mode coupling. 
These spectra may display intermittency at small scale, 
a phenomenon that would require substantially higher 
numerical resolutions to study. 

As the Reynolds numbers are increased, the amplitude 
of the nonlinear term can be expected to increase, and 
hclicity should be given by 

(Hi^)^,. = -J {'"^)^ -jiu-d^u- Vu))j_^, . 

(11) 

In the simulations presented in the following section, we 
compared the ratio of the second term to the first in 
Eq. pT|) . For simulations with 0.5 < N/f < 3 and with 
ReRo^ < 20 and ReFr"^ < 20, the amplitude of the sec- 
ond term is smaller than 10 — 20% of the first term in 
most of the runs, and increases for runs with larger values 
of N/f. 

Therefore, we can expect that for very small values 
of the control parameters, and for flows for which the 
geostrophic and hydrostatic balance holds, the latter im- 
plying dzW = 0, the helicity should remain zero. As fluc- 
tuations develop, and as weak nonlinear perturbations 
come into play, small departures from geostrophy will 
develop allowing for non-zero correlations between buoy- 
ancy and vertical velocity, as appears in Eqs. ([5]) and 
and as can also be expressed in two-scale turbulence 
closure formalisms. In that case, we can expect helicity 
to develop. Helicity should be created at large scales, 
where rotation and stratification dominate over the non- 
linear term, and may be transferred to smaller scales. 
Finally, for runs with stronger nonlinearities and large 
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values of the control parameters, we can expect devia- 
tions from the prediction in Eqs. ^ and associated 
with the extra term in Eq. ([TT|) . A detailed study of these 
deviations is left for future work. 



IV. PARAMETRIC STUDY 

We have performed nine runs on grids of 512'^ points, 
and 36 runs on 256'^ grids, up to the peak of dissipation 
and beyond, with similar (but not identical) random ini- 
tial conditions and N/f e [1/2, 16.7]. 

Figure [2] gives the temporal evolution of helicity (top 
and middle) for several runs at either fixed Fr or fixed 
Ro, and of kinetic enstrophy Zy (bottom) for several 
flows at fixed Fr; the potential enstrophy Zp shows a be- 
havior similar to Zy, with slightly smaller values. Note 
that in all quantities the oscillations are due to gravity 
waves because of the fact that our initial conditions arc 
chosen to be unbalanced, and their periods arc propor- 
tional to N. Across all runs, the maximum of Zy varies 
from 30 (for weak waves) to « 2.5, corresponding to the 
smallest Froude number considered. The time to reach 
this maximum varies from 1.5 to 3.2 tnl. The growth 
of enstrophy is typical of a turbulent flow, and is due to 
vortex stretching. The growth in the presence of waves 
is weaker, a characteristic of a weak turbulence regime. 

The overall structures in this type of flows are shown 
in Fig. [21 which displays volume rendering of buoyancy 
right after the peak of enstrophy for a run with Fr = 0.1 
and N/f = 4 (left), and for a run with Fr = 0.025 and 
N/f = 2 (right), both performed on grids of 512'^ points 
and with identical initial Reynolds numbers. The 3D ren- 
dering puts in evidence the stratification and the presence 
of large-scale layers; small-scale features with curved rib- 
bons also occur for the run with smaller stratification. 
The run with Fr = 0.1 shows strong turbulent fluctua- 
tions, whereas the run with Fr = 0.025 is smoother, with 
weaker small-scale fluctuations. 

We now examine the relation given by Eq. ([5]). In Fig.|3] 
is given the variation with the vertical index z„ (i.e., the 
vertical grid point, in runs with Up = 256) of the ratio 
r = — {H±)j^ J / {wdzO)±^f The sub-index t indicates the 
quantities were not only averaged in planes perpendicular 
to z, but also averaged in time over the peak of enstrophy. 
The ratio is shown for two runs with different values of 
N/f, namely 1.5 (top) and 16.7 (bottom); the horizontal 
line gives the prediction based on weak nonlinearities, 
i.e., N/f. One observes regular variations around the 
mean value in the vertical, so in the following we shall 
perform vertical averaging as well. These fluctuations are 
likely associated with alternating quiescent and turbulent 
patches where the advection term is strong. In spite of 
these fluctuations, the run with N/f ~ 1.5 shows good 
agreement with the prediction based on weak nonlinear- 
ities, while the run with N/f = 16.7 does not. After 
performing a detailed analysis of all the runs, it is found 
that for all runs a good agreement with Eq. ([5]) obtains 
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FIG. 4: Variations of r = - {H^)^ ^ / {wd^B) (see Eq. [8]) 
with vertical layers of index Zn, n £ [1,256] is the index 
of the vertical plane, and the data is temporally averaged 
around the peak of enstrophy. The horizontal lines give the 
geostrophic balance prediction. Both runs are performed on 
grids of 256^ points with Tie « 4 x 10^, Fr = 0.0127, and 
ReFr'^ ^ 0.672. Top: N/f = 1.5, ReRo'^ « 1.51. Bottom: 
N/f = 16.7, ReRo'^ 186.6. In the latter case, the predic- 
tion stemming from assuming weak nonlinearities no longer 
applies. 



for N/f < 3; it is also fulfilled when ReRo'^ < 20 to- 
gether with ReFr^ < 20, i.e., for strong enough waves 
and weak nonlinearities as explained in the previous sec- 
tion. Note that the importance of the buoyancy Reynolds 
number has been identified previously, e.g., in the context 
of an emphasis on the role of anisotropy and the onset 
of Kelvin-Helmholtz instabilities due to vertical shearing 

Given the measurable vertical variations observed in 
Fig. m and following the expression in Eq. ([8]) , we display 
in Fig.[5](top) a scatter plot of {Hj_{t = 0)) {H±{t)) ^ 
as a function of N/f {wdzd) for all runs with N/f < 3; 
^ =J-,t,z represents averaging on horizontal planes, for 
half an eddy turn-over time after the maximum of enstro- 
phy, and over all the vertical planes as well. This allows 
for smoothing over temporal variations due to gravity 
waves, and over the vertical inhomogeneities of the flow 
that are inherent in strongly stratified flows as discussed 
before. The symbols in Fig. [5] indicate different Froude 
numbers, and the flUcd symbols are used for runs on grids 
of 512^^ points. For the runs with N/f < 3, all points 
lie close to a straight line with slope one, showing that 
the helicity created and the source of helicity according 
to Eq. ([8]) are linearly correlated. The result presented 
in the flgure is robust for different choices of range over 
which the temporal average is performed, with windows 
of 0.5, 1 and 1.5 eddy turn-over times after the peak of 
the enstrophy. When averaging over later times, say from 
3.6 to 7.2, good agreement with the prediction of helical 
geostrophic balance also holds, in part due to the fact 
that at late times, the Reynolds number has decreased 
and waves are now more easily predominant, with smaller 
Froude and Rossby numbers. It is interesting that the 



8 



range of validity in N / f corresponds in part to the range 
identified in [53| on tlie basis of a lack of resonant inter- 
actions for these parameters. 

We also show in the shaded insets of Fig. [5] (top, mid- 
dle) the same quantity for all the runs (i.e., all values of 
iV//, corresponding to all 9 runs with 512'^ points and 
36 runs with 256^ points). For N / f > 3, creation of hc- 
licity still occurs, though not quite at the level predicted 
by Eq. ([8|). The middle graph in Fig. Ogives the same 
scatter plot but thresholded for the buoyancy and the 
inertial Reynolds numbers, TZb^tz < 20 (the inset shows 
again all points for comparison) . In both cases of thresh- 
olding, about half the points are selected (roughly, 20), 
and the points in common between the top and middle 
scatter plot are 80% that again (namely, 16, see Table 
IJ). We can interpret this fact by saying that scales (tem- 
poral and spatial) are not sufficiently separated and it 
is difficult to sort out what may be the dominant effect 
for flows to obey the relationship of Eq. ([S]) : a compara- 
ble rotation and stratification, or simply a low buoyancy 
Reynolds number. Finally, for completion, we give in 
Fig. M (bottom) {Hj_{t = 0))^ - {Hj_{t))^ as a function 
of N/ f alone. Note how most of the points pile up near 
negative values, even for large values of N/f. Since a 
growth of net helicity, and similarly of relative helicity, is 
not observed in freely decaying 3D homogeneous turbu- 
lence, with no rotation and no stratification, this further 
confirms that the production of helicity is characteristic 
of the regime under study. We finally note that it seems 
to be controlled more by the imposed stratification than 
by the Rossby number, in agreement with the fact that in 
rotating turbulence, helicity is conserved in the absence 
of dissipation. 



V. CONCLUDING REMARKS 

A parametric study of decaying rotating stratified tur- 
bulence shows that helicity is spontaneously produced at 
large scales, and that for N/f < 3 (or, ReFr^ < 20 to- 
gether with ReRo^ < 20), its value is associated with cor- 
relations between buoyancy and vertical shear, as derived 
in [2l| (for non-stratified flows see [i^l, and for the mag- 
netic case see [l^l). This creation of helicity still takes 
place for larger values of N/f, and thus confirms the pos- 
sibility, for geophysical and astrophysical flows, that the 
combination of rotation and stratiflcation creates helicity 
which in turn can be the source of large-scale magnetic 
flelds, as observed in stars and planets. 

Helicity production in rotating stratiflcd flows can also 
be related to the observation of large-scale helicity in the 
atmosphere of the Earth, although it is not occurring 
in our study through an instability involving anisotropic 
small-scale helicity as studied before in jsi], H^, but 
rather through a quasi-linearization of the large-scale dy- 
namics. Such large-scale helical flows might be relevant 
to the persistence of large-scale convective storms and to 
the onset phase of hurricanes [1, H^]. It has also been 
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FIG. 5: Helical geostrophic balance {H±{t — 0))^ ^ — 
{H±{t))^ as predicted by Eq. ((§]), temporally and vertically 
averaged for runs with N/f < 3 (top). In the shaded insets, 
the same data is given for all 45 runs of this study. Each 
symbol corresponds to a value of Fr (see labels), and the 9 
solid symbols indicate runs on grids of 512'' points at higher 
Reynolds numbers. In the middle graph, the points are se- 
lected using TZb < 20, TZn < 20 (the inset shows again all 
points for comparison). For completion the same scatter plot 
as a function of N/f only, using the same symbols, is given 
at the bottom. 
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shown that hchcal motions can be associated with the 
spiral rainbands of hurricanes when taking moisture into 
account in the dynamical equations [3ll |. 

The observed saturation in the level of helicity for 
larger values of N/ f and for sufficiently strong stratifi- 
cation can likely be understood in terms of the presence 
of vertically-sheared horizontal flows in that regime (see, 
G-g-, [43,7 531), a tendency that persists in the absence 
of rotation [5^. The generation of helicity requires an 
interplay between stratification and rotation, and when 
stratification dominates, vertical and horizontal motions 
are less correlated. For large values oiReFr^ and ReRo^, 
the deviations from the prediction assuming weak non- 
linearities can also be associated with an increase of the 
amplitude of the nonlinear term in the momentum equa- 
tion, as verified in the simulations by direct estimation of 
the amplitude of the different terms in the equation for 
the helicity. 

The fact that rotating stratified fiows can sponta- 
neously produce large-scale helicity opens new lines of 
research and questions. For example, is there a detailed 
role to be played by potential vorticity conservation on 
the emergence of helicity? How would the inclusion of ei- 
ther shear, radiation, moisture, or some general forcing in 
Eq. ([3|) modify these results? And finally, how would tur- 
bulence affect significantly the creation of helicity, as the 
Reynolds numbers are further increased? Indeed, mixing 
is thought to have two transitions in terms of TZ r- i n the 
presence of an imposed shear, it was shown in [55| (see 
also [iBl) that below 7, molecular diffusion is observed. 



with basically no turbulence; the intermediate regime 
7 < TZb < 100 follows a linear Osborn diffusion law [56l |. 
and above that value, a new regime is reached with dif- 
fusivity scaling as Ti-^J^ ■ The latter regime is of course 
what matters for geophysical fiows with TZb ~ 10*, such 
as in the meridional overturning circulation, central to 
climate dynamics, but is it affected by the production of 
helicity? 

We thus plan to pursue this study concerning the role 
of helicity in rotating stratified turbulence in order to 
help decipher the different mechanisms at play, a study 
that will eventually lead to better sub-grid scale models 
of such flows that arc needed to obtain a more accurate 
representation of enhanced diffusivities in weather and 
climate models. 
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